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Abstract. We study preconditioning techniques for discontinuous Galerkin discretizations 
of isotropic linear elasticity problems in primal (displacement) formulation. We propose 
, subspace correction methods based on a splitting of the vector valued piecewise linear dis- 

' continuous finite element space, that are optimal with respect to the mesh size and the 

Lame parameters. The pure displacement, the mixed and the traction free problems are 
discussed in detail. We present a convergence analysis of the proposed prcconditioners 
and include numerical examples that validate the theory and assess the performance of the 



O 

' prcconditioners. 
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1. Introduction 



< 

^ . The finite element approximation of the equations of isotropic hnear elasticity may be 

, accomphshed in various ways. The most straightforward approach is to use the primal for- 
mulation and conforming finite elements. It is well known that such a method, in general, 
does not provide approximation to the displacement field when the material is nearly incom- 
^ . pressible (the Poisson ratio is close to 1/2). This phenomenon is called volume locking. To 
[ alleviate locking, several approaches exist. Among the possible solutions, we mention the use 
^ I of mixed methods, reduced integration techniques, stabilization techniques, nonconforming 
■ methods, and the use of discontinuous Galerkin methods. We refer to [11, 14] for further 
discussions on such difficulties and their remedies. In this work we focus on the Symmetric 
Interior Penalty discontinuous Galerkin (SIPG) methods introduced in [14, 15, 19, 20] for 
the approximation of isotropic linear elasticity. We have chosen to work with these DG dis- 
cretizations, since we have in mind a method that is simple but still applicable to different 
types of boundary conditions. In fact, unlike classical low order non-conforming methods 
' ^^^^ [H])' Interior Penalty (IP) stabilization methods introduced in [14, 15] can be shown 

^ ' to be stable in the case of essential (Dirichlet or pure displacement) boundary conditions, 
or natural (Neumann type, or traction free) boundary conditions. As a consequence, these 
IP methods provide a robust approximation to the displacement field and avoid the volume 
locking regardless the boundary conditions of the problem. 

For the design of the preconditioners we follow the ideas introduced in [4] for second 
order elliptic problems. However, such extensions are not straightforward, since we aim at 
constructing preconditioners that work well for three different types of boundary conditions: 
essential, natural and mixed boundary conditions, used in linear elasticity. This complicates 
the matters quite a bit. We consider a splitting of the vector valued, piecewise linear, 
discontinuous finite element space, into two subspaces: the vector valued Crouzeix-Raviart 
space and a space complementary to it which consists of functions whose averages are 
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orthogonal to the constants on every edge/face of the partition. This space decomposition 
is direct and the spaces are orthogonal with respect to a bilinear form obtained via using 
"reduced integration" to calculate the contributions of the penalty terms in SIPG. 

In the pure displacement case (essential boundary conditions), the restriction of the bi- 
linear form based on reduced integration is coercive on the Crouzeix-Raviart space and is 
spectrally equivalent to the SIPG bilinear form. The space decomposition mentioned above is 
then orthogonal in this reduced integration bilinear form. Thus, in case of essential boundary 
conditions we have a natural block diagonal preconditioner for the linear elasticity problem: 
(1) a solution of a problem arising from discretization by nonconforming Crouzeix-Raviart 
elements; (2) solution of a well-conditioned problem on the complementary space. 

For traction free problems or problems with Dirichlet conditions only on part of the bound- 
ary, the situation is quite different. On one hand the reduced integration bilinear form when 
restricted to the Crouzeix-Raviart space has a null space whose dimension depends on the 
size of the problem (see [11]). On the other hand in the full SIPG bilinear form (without 
reduced integration) the space splitting discussed above is no longer orthogonal. Our ap- 
proach in resolving these issues is based on a delicate estimate given in §3.1 which shows a 
uniform bound on the angle between the Crouzeix-Raviart and its complementary space in 
the SIPG bilinear form for all types of boundary conditions. Once such a bound is available 
we show that a uniform block diagonal preconditioner can be constructed. 

The rest of the paper is organized as follows. We present the linear elasticity problem, the 
basic notation and discuss the DG discretizations considered in §2. Next, in §3 we introduce 
the splitting of the vector valued piecewise linear DG space and discuss some properties of 
the related subspaces. In section §4, we introduce the subspace correction methods, and we 
prove that they give rise to a uniform preconditioner for the symmetric IP method. The last 
section §5 contains several numerical tests that support the theoretical results. 

2. Interior Penalty Discontinuous Galerkin methods for linear elasticity 

equations 

In this section, we introduce the linear elasticity problem together with the basic notation 
and the derivation of the Interior Penalty (IP) methods and we discuss the stability of these 
methods. 

2.1. Linear Elasticity: Problem formulation and notation. Let C IR'^, d = 2,3, be 
a polygon or polyhedron (not necessarily convex) and let n be a vector field in IR'^, defined on 
Q such that u G [H^{Q)Y. The elasticity tensor, which we denote by C, is a linear operator, 
i.e., C : IRgym ^ ^sym^ acting on a symmetric matrix A G IRfym? in the following way: 

C A = 2/iA + Atrace(A)J, 

where /i and A are the Lame parameters and satisfy < /ii < /i < /i2 and < A < oo. 
In terms of the modulus of elasticity (Young's modulus), CB, and Poisson's ratio, u, the 
Lame parameters can be rewritten in the case of plane strain as: fi = €/{2{l + u)) and A = 
i/€/((H-z/)(l — 2z/). The material tends to the incompressible limit (becomes incompressible) 
when the Lame parameter A — )■ oo or equivalently when the Poisson's ratio z/ — )■ 1/2. 

One can show that the linear operator C is selfadjoint and has two eigenvalues: (1) a 
simple eigenvalue equal to (2/i -|- d\) corresponding to the identity matrix; (2) an eigenvalue 
equal to 2/^, corresponding to the ^Mlll _ i dimensional space of traceless, symmetric, real 
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matrices. Thus for (i = 2, 3, we always have that 

(2.1) 2/i(A : A) < {CA : A) < {2fi + dX){A : A), 

where (■ : ■) denotes the Frobenius scalar product of two tensors in M!^^'^. We also denote 
by (■, ■) the Euclidean scalar product of two vectors in IR"', i.e., 

d d d 

{v, w) = ^ VkWk, (v : lu) = ^ ^ VjkWjk- 

k=l j=l k=l 

The corresponding inner products in [L^(f2)]'^ and [L'^{Q)Y^'^ are denoted by 



{v,w) = / {v,w), {v : w) = {v : w). 
Jn Jn 

We write dQ = UT^i with Fjv and referring respectively to the subsets of the dQ 
where Neumann and Dirichlet boundary conditions are imposed. 

Let s{u) = |(V'U + (Vn)^) be the symmetric part of the gradient of a vector valued 
function u. The elasticity problem in primal formulation then is: Find u G [Hy'^"{Q)Y, 
a > 0, which is the unique minimizer of the energy functional J^{u), given by 

(2.2) J{u) := ^iCe{u) : £(n)) - (/, u) - (g^, v)r^ 

Here / G [//^(fi)]*^ is a given volume force and qn G [H^^'^{T]\f)Y is a given surface force 
acting on F^r C dQ. The Euler-Lagrange equations corresponding to the minimization prob- 
lem (2.2) give the following well known system of linear PDEs for the unknown displacement 
field u: 

— div(C£(n)) = /, on 

(2.3) {Ce{u))n = qn, on Fat, 

u = 0, on Vr,. 

In the above equations, n is the outward unit normal vector to dVt. The solution u vanishes 
on a closed part of the boundary Td (Dirichlet boundary) and the normal stresses are 
prescribed on Fa? (Neumann part of the boundary). In the traction free case (F^r = dVL), 
the existence of a unique solution to (2.3) is guaranteed if the data satisfy the following 
compatibility condition: 



f vdx+j gN-vds = G RM(fi), 
n Jan 

where RM(f2) is the space of rigid motions, defined by: 

(2.4) RM{n) := {v = a + bx : aeR'^ be so{d) } 

where x is the position vector function in Q and 5o{d) is the Lie algebra of skew-symmetric 
dxd matrices. In this case, the uniqueness of solution is guaranteed up to a rigid motion (and 
is unique, if we require that the solution is orthogonal to any element from RM(f2)). In the 
case of F^) 7^ and closed with respect to dQ no extra conditions are required to guarantee 
uniqueness. By considering the variational formulation of (2.3), the issue of solvability and 
uniqueness of the problem reduces to show coercivity of the associated bilinear form. As it 



4 



BLANCA AYUSO DE BIOS, IVAN GEORGIEV, JOHANNES KRAUS, AND LUDMIL ZIKATANOV 



is well known, for linear elasticity, this hinges on the classical Korn's inequality [10] which 
guarantees the existence of a generic positive constant Cn > such that: 

(2.5) \\Vv\\l^ < Cn {Mv)\\l^ + \\v\\l^) , V^; G [H^n)^ • 

The second term on the right hand side can be omitted as follows from the Poincare or 
Poincare-Friedrich's inequality, obtaining thus first Korn's inequality for v G [i/g p^(r2)]'^ 
and second Korn's inequality for v G [H^{Q)Y /IllSA{Q). 

2.2. Interior penalty methods: Preliminaries and notation. We now introduce the 
basic notations and tools needed for the derivation of the DG methods. 
Domain partitioning. Let Th be a shape-regular of partition of Q into d- dimensional 
simplices T (triangles if = 2 and tetrahedrons if d = 3). We denote by hx the diameter 
of T and we set h = maxx^Th ^t- We also assume that Th is conforming in the sense that it 
does not contain hanging nodes. A face (shared by two neighboring elements or being part 
of the boundary) is denoted by E. Clearly, such a face is a {d — 1) dimensional simplex, that 
is, a line segment in two dimensions and a triangle in three dimensions. We denote the set 
of all faces by Sh, and the collection of all interior faces and boundary faces by and £f , 
respectively. Further, the set of Dirichlet faces is denoted by Sj;', and the set of Neumann 
faces by Sj^ . We thus have. 

Trace operators (average and jump) on E G Sh- To define the average and jump trace 
operators for an interior face E G Sj^, and any T G Th, such that E G dT we set Ue.t to 
be the unit outward (with respect to T) normal vector to E. With every face E & Sj^ we 
also associate a unit vector Ue which is orthogonal to the (d — 1) dimensional affine variety 
(line in 2D and plane in 3D) containing the face. For the boundary faces, we always set 
Ue = TiE^Ti where T is the unique element for which we have E C ST. In our setting, for 
the interior faces, the particular direction of is not important, although it is important 
that this direction is fixed. For every face E G Eh-, we define T^{E) and T~{E) as follows: 

T+{E) := {T eTh such that E C dT, and {ue, Uet) > 0}, 

(2.6) 

T-{E) := {T eTh such that E C dT, and {ue, Ue^t) < 0}. 

It is immediate to see that both sets defined above contain no more than one element, that 
is: for every face we have exactly one T~^{E) and for the interior faces we also have exactly 
one T~{E). For the boundary faces we only have T~^{E). In the following, we write 
instead of T^{E), when this does not cause confusion and ambiguity. 

For a given function w G [L^(r2)]'^ the average and jump trace operators for a fixed E ^ 
are as follows: 

(2.7) {{w} := (^^^^-^) , M := (^+ " w'), 

where and w~ denote respectively, the traces of w onto E taken from within the interior 
of T~^ and T~. On boundary faces E G Sj^, we set = w and {wj = w. We remark 

that our notation differs from the one used in [1], [3], [2] (which is considered a classical 
one for the IP methods). We have chosen a notation that is consistent with the one used in 
[15], where the IP method we consider was introduced for the pure displacement problem. 
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In addition, it seems that such a choice leads to a shorter and simpler description of the 
preconditioners we propose here. 

Finite Element Spaces. The piecewise linear DG space is defined by 

V^^ ■= {u G L'^in) such that G ¥\T), VT G % }, 

where P^(T) is the space of linear polynomials on T. The corresponding space of vector 
valued functions is defined as 

For a given face E, we denote by : L'^{E) h-t- P°(£') the L^-projection onto the constant 
(vector valued or scalar valued) functions on E defined by 

(2.8) Vlw = I w for all w G L^{E), 

\E\ Je 

(2.9) V%w = [ w for all we [L^E]]'^. 

\E\ Je 

Observe that for w G the mid-point integration rule implies that V'^w = winiE) for 
all E G £h, with uie denoting the barycenter of the edge or face E. 

The classical Crouzeix-Raviart finite element space can be defined as a subspace of V^^ , 
as follows: 

(2.10) l^^^ ={vE V^"^ : Vl[v] =0, \/Ee 81] . 
The corresponding space of vector valued functions is 

(2.11) V^^ := [V^^Y 

2.3. Weighted residual derivation of the IP methods. In [15] the authors introduced 
a symmetric interior penalty method for the problem of linear elasticity (2.3) in the pure 
displacement case (i.e, = dVL, = 0). We define the function space 

[H^iTh)]'' ={ue [L^n)]'' such that G [H^T)]'', VT G } . 

For any pair of vector fields (or tensors) v and w, we denote 

{v,w)r,^ = ^ / {v,w). 

Ten 

For scalar and vector valued functions we also use the notation 

(2.12) [v,w)£ = Y^ / vw, and {v,w)£ = Y^ / {v,w) . 

E&s"^^ Eee"^^ 

We now derive, using the weighted residual framework [8], the IP methods for the more 
general case of mixed boundary conditions. To present a short derivation of the methods, 
we assume u G [H'^{fl)Y. Such assumption is not required for the methods to work. We 
present the derivation under such assumption in order to avoid unnecessary details which 
would shift the focus of our presentation on preconditioners. 
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-div(C£(n)) 


= f 


on T eTh 


l{Ce{u))n\E 


= 


onE e£l 


Me 


= 




Me 


= 


on E eSj^ 


J{Ce{u))n-gN]E 


= 


on EeS^ 



By assuming that the solution of (2.3) is a priori discontinuous, u G [H'^{Th)Yi we may 
rewrite the continuous problem (2.3) as follows: Find u G [H'^{Th)Y such that 



(2.13) 



where we recall that Ce{u) = 2fie{u) + A trace(£(w))/. Following [8], we next introduce a 
variational formulation of (2.13) by considering the following five operators 

and weighting each equation in (2.13) appropriately. This then amounts to considering the 
following problem: Find u G [H'^{Th)Y such that for all v G [H'^{Th)Y 

-diy{Ceiu)) - f,Boiv))r, + (I(Cs(u))n], S2(t;)),o + (M, i3i(i;)),o 



Bo 
Bi 

B2 



(2.14) 



{lulBfiv)),n + iliCeiu))n-g4,BUv)),. = 0. 



Different choices of the operators Bq, Bi, B2, Bf and B2 above give rise to different variational 
formulations and, consequently to different DG methods. We refer to [8, Theorem 6] for 
sufficient conditions on the operators Bq, Bi, B2, Bf and to guarantee^ the uniqueness of 
the solution of (2.14). 

To derive the IP method of interest, we take v piecewise smooth and we set Bq{v) = v, 
B2{v) = ivj and (v) = v in (2.14), to obtain that 



(2.15) 



{-div{Ce{u)),v)r, + {liCeiu))n],iv}}),o^,o + {lulB^iv)),o^,o 
= if,v)r^ + i9N,v)sN . 
Defining 

(2.16) nv) = {f,v)r, + {lgll3f{v)),n + ig^,v),., 

and integrating by parts the first term on the left side of (2.15) then leads to 

(2.17) {Ce{u) : e{v))r, - iUCe{u))n}}, M)eoue- + iM,Mv)),oue- = H^)- 
For a fixed edge E E Sj^U Ej^ the operator Bi{v) is defined by 

(2.18) B,{v) := -{{(Cs(^))n}} + a,P,V%[v\ + ai/3iH, 

where, following [15], the parameters (3o and (5\ are chosen depending on the Lame constants 
A and /i: 



(2.19) 



/3o := d\ + 2/i, 



/3i := 2/i . 



""^We note that in [8] the focus is on the scalar Laplace equation. The arguments for the elasticity problem, 
are basically the same. 
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The remaining two parameters, and ai, are still at our disposal to ensure (later on) 
stability and to avoid locking of the resulting method. 
We define 



(2.20) 



and set 



a,,o(M,H):=ao/3o /(^^'M'^M^ 



a,(M, M) = a,-o(M, M) + a,, i(M, H). 
Then, the weak formulation of Problem (2.13) reads: Find u G [H'^{Th)Y such that 

(2.21) A{u, w) = T{w), ywe [H\rh)Y. 
The bilinear form .4(-, ■) is given by 

(2.22) A{u, w) = A(ti, + a,, i(|n], H), 
where 

^2 23) A(n,ti^) = {Ce{u):e{w))r,-{l{Ce{u))n^,lw]),o^sf; 

-(M, §(Cs(^))n}}),.u.f + a,,o(M, M). 
It is straightforward to see that 

^2 24) ^(^'^) = {Ce{u):e{w))r,-{i{Ce{u))n^,lwl),o^,u 

To obtain the discrete formulation, we replace the function space [H'^{l'h)Y (2.21) by 
and we get the IP-1 approximation to the problem: Find G such that: 

(2.25) A{uh, w) = J^{w), \/ we 

We could also consider the approximation given by the IP-0 method: Find Uh G such 
that: 

(2.26) Ao{uh, w) = T{w), V G V^^. 

As we see next, the IP-0 method provides a robust approximation to the problem (2.3) 
in the pure displacement problem To = dQ. As we mentioned earlier, for other types of 
boundary conditions such equivalence in general does not hold. 

Remark 2.1. Although we do not consider non- symmetric IP methods in this paper, let us 
remark that non- symmetric versions can be easily incorporated in the definition of Bi{v). 
For example, by setting: 

B,iv) := 9{{iCeiv))n}} + aoPoVU^j + «i/3iH, 

we obtain a non-symmetric bilinear form for the values 6 = or 6 = 1. Such values of 6 
correspond to the Incomplete Interior Penalty (IIPG, 6 = 0) and Non-symmetric Interior 
Penalty (NIPG, 9 = 1) discretizations, respectively. 
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2.4. Stability Analysis. We close this section presenting the stability and continuity results 
pertinent to our work. We start by introducing some norm notation. For v G [H^(Th)Y "^^ 
define the semi-norms 



iVt^llk 



(2.27) ^^^^ 



|2 



and norms: 



(2.28) 11^11^= ||C^/'s(^)||^,-^+/3o|P^MII^ + AINI^+ E hE\\C'/Mv)-n\ 



0,E ■ 



Ee£°^U£^ 



For V G V^'^ we define the norms 

(2.29) IMlLGo = l|C^/'^WIIk+/3ol^M«ll* 
and 

(2.30) \\v\\lG = Mlco + mA\l- 

Notice that for v G V^*^ the norms (2.28) and (2.30) are equivalent. We finally introduce 
the norm: 

(2.31) M'miT,) = l|V«||^,r, + /3o|P^MII^ + +/3ilNII^ . 

Notice that continuity of the IP-1 and IP-0 bilinear forms with respect to the norm (2.28) 
follows easily from Cauchy-Schwarz inequality together with the bound on the maximum 
eigenvalue of C, i.e., for all u G [H'^{Th)Y ^^"^ ^ V^'^ we have 

(S(C£(n)n)S,H),.,,. = {i{Ce{u)n)}},V%[v\ 

OlQ ' h h n 

The equivalence of the norms (2.28) and (2.30) for any v G V^'^ guarantees therefore the 
continuity of the IP-1 bilinear form with respect to the norm defined in (2.30) for finite 
element functions. 

The solvability of the discrete methods (2.25) and (2.26) is guaranteed if and only if, a 
discrete version of the Korn's inequality holds on V^'~^ . In [7] the following discrete Korn 
inequality is shown for [if^(7h)] ''-vector fields: 

(2.32) W^AIt, < C {\\e{v)\\lr, + \'niM\l + || V x ^H^^^-J 

where tti : [L?'{£h)Y — ^ "^^{^h) is the L^-orthogonal projection onto the space of piecewise 
linear vector valued functions on £h (or a subset of it). 

Coercivity of the IP-1 bilinear form with respect to the norm (2.30) can be easily shown 
by taking u = w = v in (2.24): 

A{v,v) = {Cs{v) : e{v))r,+aoMh-E'^'vUv}\\leour,+ai^^^^^^^ 
-2i{{iCeiv)n)},lv]),.^,n . 
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Using Cauchy-Schwarz, trace and inverse inequalities together with the arithmetic-geometric 
inequahty and the bound on the maximum eigenvalue of C it follows that 



{i{Ce{v)n)}},M)eo^e- = (S(Ce(t^)n)}}, Pj^H),. 



D 



+ C'mi)) 2 , C^O/^O II ,-1/2^0 ir 11112 

«oPo ^ 



(2.33) < + + ^l|/^^^/XM||g,,.,,,. 



Hence, we finally have 



A{V^V) > (l_ 2^(1 + an.) )||^l/2^(^)||2^^^ ^^^^11^-1/2^^,112 



"0 



2 

and therefore by taking ao = max(l,4Cf(l + Cinv)) (sufficiently large) we ensure the coer- 
civity of A{-, ■) with respect to the || ■ Hbg-hoi'iii with constant independent of h, fi, and A. 
Using now (2.32) (since the norm (2.30) contains the full jump) we conclude that A{-, ■) is 
coercive with respect to the || ■ ||/^i(7yj-norm (2.31). Therefore the IP-1 method defined by 
(2.24) provides a robust approximation to (2.3) and does not lock as A — oo. 

As we mentioned earlier, in the pure displacement case {Td = dQ, Ftv = 0) the bilinear 
form Ao{--) defined in (2.23) is coercive. Indeed we may use the identity (which holds for 
C^{n) functions): 

(2.34) div£(v) = ^ (divVv + Vdivv) 

and rewrite the volume term in (2.23) (also in (2.24)) as follows: 

{Ce{u) : e{w))r, = E / (^^(^) ■ ^(^)) 

Ten -^"^ 



E j (2/i(V'U : V-y) + (/i + A)(divw,divt;)) . 



Then, from the discrete Poincare inequality [12, 6], the resulting modified bilinear form for 
^o(') ■) is now coercive in with respect to the || • ||j|^i(7^j norm, with coercivity constant 
independent of h and A; 

(2.35) A^{v,v)>C\\v\\]j,^j-^^ V^gV^^. 

Therefore, the discrete problem (2.26) is well posed and the IP-0 method is stable and robust 
(locking free in the limit A — >■ oo). Notice that in (2.35) we are using the || ■ ||j|/i(7-^)-norm 
which includes not only the norm |Pg[t)]|^,, but also the norm Hi']]*. This is a consequence 
of the vector valued counterpart of [4, Lemma 2.3]. The stability property given in (2.35) 
implies that the IP-0 and IP-1 methods are spectrally equivalent for the pure displacement 
problem. These observations are summarized in the next Lemma: 

Lemma 2.2. Let ^(•, •) and Aq{-, ■) he the bilinear forms of the IP-1 and IP-0 methods 
for the linear elasticity problem, defined in (2.24) and (2.23), respectively. For the pure 
displacement problem Fx) = dVL, r^r = 0, there exist a constant c > that depends only on 
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the geometry of the domain f2 but is independent of the mesh size and the Lame parameters 
fi and A such that 

(2.36) Aoi^, < < cAoiv, v) Wve V^^ . 

The above lemma guarantees that for the pure displacement problem, constructing a 
uniform preconditioner for the IP-1 is equivalent to constructing a uniform preconditioner 
for the IP-0 method (see [4]). For linear elasticity equations, unlike for scalar equations, 
this can be done only when Tu = dQ. 

For a detailed derivation and error estimates, we refer to [15, Theorem 2.5]. 

3. Space decomposition 

We present now a decomposition of the DG space of piecewise linear vector valued func- 
tions that plays a key role in the construction of iterative solvers. This decomposition was 
introduced in [4] for scalar functions and also in [9] in a different context. Its extension to 
vector valued functions is more or less straightforward. We omit those proofs which are just 
an easy modification of the corresponding proofs in the scalar case. However, we review 
the main ingredients and ideas behind such proofs, since they play an important role in the 
analysis of the preconditioner given later on. In the last part of the section we give some 
properties of the spaces entering in the and prove a result that is essential for showing that 
the proposed preconditioner is uniform. 

Following [4] we introduce the space complementary to V'^^ in V^*^, 

(3.1) Z = {ze V^^ and = 0, for all E e Sj^} . 
The corresponding space of vector valued functions is 

(3.2) Z = \Zf. 

To describe the basis functions associated with the spaces (2.11) and (3.2), let '~Pe,t denote 
the scalar basis function on T, dual to the degree of freedom at the mass center of the face 
i?, and extended by zero outside T. For E e (9T, E' G 9T, the function if>E,T satisfies 



1 \i E = E', 
otherwise, 



and also we have 
For all u G we then have 

(3.3) U{x) = UT{mE)^E,Ti^) = XI U^{^E)'fE{^) + ^ u~ {mE)^ 



E\X)^ 



where in the last identity we have just changed the order of summation and used the short 
hand notation ^%{x) := (pE,T±ix) together with 

u^{mE) := UT±{mE) = ^ / UT±ds, \/ E e S^, ■ E = dT+ f] dT' , 



E 



E 



1 



u{mE) ■= UTirriE) = 7-^^ / uxds, \/ E & 8^, such that E = dT fl dVt.. 
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Recalling now the definitions of T^{E) and T~{E) given in (2.6) we set 



(3.4) 


T E 


= ^E,T+(E), 


- 09 7? T'- 


(E)} 




and 














Ve- 


^E,T+{E) — 


^E,T- 




WEeS^. 


(3.5) 


2 

= ^E,T+{E), 







Some clarification is needed here. Note that from the definition of (Pe,t+[e) cind (Pe,t-{e) 





Figure 3.1. Basis functions associated with the face E: (left) and (p^^ 
(right). 

for an interior edge E e 8^, it does not follow that their sum is even defined on the edge E, 
since it is just a sum of two functions from L'^{fl). However, the sum {^e,t+{e) + Ve,t-{e)) 
has a representative, which is continuous across E and this representative is denoted here 
with 93^^, see Figure 3.1. 

Clearly, {'{>'e^} E&efvje^ linearly independent, and {'^Pe} E&£°u£f linearly indepen- 
dent. A simple argument then shows that 

= span {W'i'^ek}t=i}j,^^o^^N , 2 = span {{tjjEekjLi} Eee^usj^ ■ 

Here e^, k = 1, . . . , d is the k-th canonical basis vector in IR"^. Hence by performing a change 
of basis in (3.3), we have obtained a "natural" splitting of 

yDG ^ yCR ^ 2 

and the set 

(3-6) {'»/'l;}sGf°Uff U {v9i'^}£;g£;ou£:Ar, 

provides a natural basis for the DG finite element space. This is summarized in the next 
proposition. 

Proposition 3.1. For any u G V^"^ there exist unique v G V'"^ and a unique z E 2 such 
that 

(3.7) . ^ . + . an, " ^ ^^"^'^ " 
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The proof of the above result follows by arguing as for the scalar case in [4, Proposition 
3.1], but proceeding componentwise. The next Lemma shows that the splitting we have 
proposed is orthogonal with respect to the inner product defined by ^o('; ■)■ 

Lemma 3.2. The splitting (3.7) V^^ = V'"^ (B Z is Aq- orthogonal. That is 

(3.8) Aq{v,z)=Aq{z,v) = {} 'iveV^^, 'izeZ. 

The proof follows straightforwardly by using the weighted residual formulation (2.15)- 
(2.23) and the definition of the spaces V'~^^ and Z. 

3.1. Some properties of the space Z. We now present some properties of the functions 
in the space Z. We start with a simple observation. From the definition of the spaces 
and Z it is easy to see that 



E l|V^llo,T = (W,§V;z}})^ou^z,. 



Applying the Schwarz inequality, one then gets the following estimate 



Y,\\yz\\l^<c\\h-"'PlM 



which is a straightforward way to see that the restriction of the IP-1 and IP-O-bilinear 
forms (even for ^ = 0, 1 as in Remark 2.1) to the space Z are coercive in the || ■ ||//i(-7-^)-norm 
(2.31) (regardless whether the boundary conditions are Dirichlet, Neumann or mixed type). 
Therefore the resulting stiffness matrices are positive definite. 

The next result provides bounds on the eigenvalues of ^o('5 ■) and A{-,-), when restricted 
to Z. 

Lemma 3.3. Let Z he the space defined in (3.2). Then for all z G Z, the following estimates 
hold 

(3.9) h-^z\\l<A^{z,z)<h~^z\\l, 



and also, 



10 ; 



(3.10) [(ao)/3o + ai(3^]h~'\\z\\l < A{z, z) < [ao/^o + a,{3,]h- 
where (3q and (3i are as defined in (2.19). 

Proof. Arguing as in [4, Lemma 5.3] (but now componentwise for vector valued functions) 
one can show that (due the special structure of the space Z). 

(3.11) h-'\\z\\l< ^^'r^Wllo,i.</^-'ll^llo- 

From the coercivity of it follows then 

aoPoh-'\\z\\l<aof3o H^M-^l llo,s < A(2, 2;) . 
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Similarly, the LF'{£h) stability of the projection together with the coercivity of A gives 
{a,P, + a,P,)h-^z\\l<a,P, Yl ^^1^^ W llo,i. + "i/^i E ^^1^^WIIo,e 

<Aiz,z), 

and so, the lower bounds in (3.9) and (3.10) follow. We next show the upper bound in (3.9), 
and the upper bound in (3.10) is obtained in an analogous fashion. Using (2.33) together 
with (2.1) we get 

Ao{z,z)<ao^o Y ^^'II^^WIIo,i.+ l|C^/'^(-^)|lk 

< /3o («o||/^^'/X Wll^o.r, +C||e(z)||^,^^) . 

Hence, the upper bound in (3.9) follows in a straightforward fashion using the trace and 
inverse inequalities together with the obvious inequality ||£:(2)||o,rh — ll^^llo,77i- D 

We close this section with establishing a uniform bound on the angle between V'"^ and 
Z in the inner product given by the bilinear form A{-, ■). The estimate is given in Propo- 
sition 3.4. It plays a crucial role in bounding the condition number of the preconditioned 
system. 

We remind that E ^ £h denotes a. {d — l)-dimensional simplex (a face), which is either 
the intersection of two (i-dimensional simplices T G 7h or an intersection of a rf-dimensional 
simplex T G 7^ and the complement of il, i.e., E = T H (IR'^ \ ^)- In the former case, the 
face E is called an interior face and in the latter it is called a boundary face. 

The proof of Proposition 3.4 requires arguments involving the incidence relations between 
simplices T eTh and faces E E £hi and estimates on the cardinality of these incidence sets. 
For the readers' convenience, we provide a list of such estimates below. 

• We define Mq{E) to be the set of dimensional T eTh simplices that contain E: 

Mo{E) ■= {T G Th: such that E eT} 

By definition, for the cardinality of this set we have \Mq{E)\ = 2 for the interior faces 
and \J\fo{E) \ = 1 for the boundary faces. 

• We define the set of neighbor (or neighboring) faces N'i{E) to be the set of faces 
which share an element with E: 

Mi{E) := {E' E Sh, such that AToiE) n A/'o(S') ^ 0} 

From Proposition A.l (see Appendix A) we have that \N'i{E)\ < (2d + 1). 

• Next, we define N'2{E) to be the set of faces which share at least one neighboring 
face with E: 

J\f2{E) := {E' E Eh. such that Mx{E) n Mx{E') ^ 0} 
From Proposition A.l we have the estimate \M2{E)\ < {2d + 1)^. 
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• For the basis functions {"ipE:} EG£°ue^ have the following relations: 



(3.12) 




and [■?/'|;](a:) = 1, for all x e E, 



(3.13) 1 1^11 (a;) I < 1, for all x G E', and all E' e X2iE). 

The above relations all follow from the definition of iPe{^) ^^"^ ^^e fact that lip^} is 
linear function on every face in £h, and therefore /^[V^l;/] = |-E||^/'|;/](m£;). 
• Finally, for E E £h, E' E Sh, and E" e £h it is straightforward to see that we have: 

(3.14) If E^^^^{E')n^^^{E") then /" = 0. 

Je 

An easy consequence from the definitions then is the following: 

(3.15) If E'iM2{E") then /" = 0, for all E e Eh- 

Je 

We finally give Proposition 3.4. To avoid unnecessary complications with the notation, we 
state and prove the result for scalar valued functions. The proof for vector valued functions 
is easy to obtain, and with the same constant, by just applying the scalar valued result 
component- wise. 

Proposition 3.4. The following inequality holds for z E Z: 

(3.16) \\hE'\M-nizmiE<{-^--) E ii^^'^'wiio,i., 

with a constant p > 1 which depends on the shape regularity of the mesh. 
Proof. Since is the orthogonal projection on the constants, we have that 

(3.17) \\h-E'^\i4-'Pum\i,E = wh-E'^'imiE - \\hE'^"vuzm,E. 

Let z E Z, i.e., z = YliE'&£°vje'^ '^e''>Pe'- From (3.12) we have that 'P^'IV^I;/] = 5ee'^ and 
hence, we may conclude that 

\\h-E"V%lz\\\l^ = E^^^'^^^^^' 

= E ^EEZI= {WZ^Z). 
E€£°yj£^ 

Here we have denoted by D : IR'^*"' i— )■ IR'^'^' a diagonal matrix with non-zero elements 
"^EE '■= and by 2 G H'^'"' the vector of coefficients z = {zE}Ee£h expansion of 

z E Z via the basis {ipE}Ee£h- 
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-V2n 



\0,E 



E'e£°ue^ 



0,E 



Eee^yjs^ 



E /'"i- E E 



ze'Ze"IWe'IWe" 



E'es°ue^ E"e£°u£^ 
ze'Ze" I E / ^e'I 



E E 

E E Ze'Ze"^E'E" = {^z,z). 



'^E'\We" 



E'es-ue^ E"ee-u£^ 



Here, § : M''^''' i— >■ K'^''' denotes the symmetric real matrix with elements 
(3.18) Se'e" = E 



E&£°U£^ 



hE^hpE'n'yE"! 



E 

EeAfi{E')fWiiE") 



hE^bpE'!ii'yE"l 



In the last identity above, we have used (3.14). Note that according to (3.15), if E' ^ N'2{E") 
then 'E'E'E" = 0. Thus, 

(§5,5)= ^ ^ ze'Ze'Se'e"- 

E'&£°\j£f^ E"&Af2{E) 



From this identity and (3.12) and (3.13), we obtain that 

\E 



\EE'E"\<\Ari{E')nAri{E")\ max —- 

E&Mi{E')nMi{E") He 



E\ 



Introducing 



p = sup -— — rr- 



< (2d+ 1) max — - 

EeMi{E')fWiiE") He 



-V2Sin)-l/2,7, ,7, 



sup 



w), w) 



{w, w) 



we obtain that 

(3.19) (§5, Z) = (D-l/2gp-l/2pl/2~^ < ~y 

This inequality can be rewritten as ^(S5, 5) < (©5, 5) and hence 

(§5, 5) - (D5, 5) < (§5, 5) - i(§5, z) = (1- -)(S5, 5). 

P P 

Note that (3.17) implies that 

(3.20) (§5,5) = (D5,5) + J2 WhE^^'iM-'PEM 



\o,Ey 



E&£h 



and thus (§5, 5) > (D5, 5). This shows that p > 1 in (3.19). 

It remains to show that p can be bounded by quantities depending only on the shape 
regularity of the mesh. Again, by (3.15) we have that: if E' ^ M2{E") then E>e'e" = 0. 



16 BLANCA AYUSO DE BIOS, IVAN GEORGIEV, JOHANNES KRAUS, AND LUDMIL ZIKATANOV 



Hence: 



< ||D~^/2§D~^/2||£oo < max V 

E'eM2{E") 

\^E'E 



iE"E' 



E"&S°VJ£S ^,.t^^,n ^^E'E'^E"E" 



< max 



E"&£°ue^ 



\M2{E")\ max 



E'eM2{E") a/I 



'E'E'^E"E" 



< {2d +1)^ max max max 



\E\ / hpih 



E" 



E"e£°ue^ E'€j\f2{E") E£j\fi{E')fWi{E") He y l-E'll-E"! 

The quantity on the right side of this estimate only depends on the shape regularity of the 
mesh and the proof is complete. □ 

Remark 3.5. We remark that the constants in Proposition 3.4 can be sharpened, at the price 
of further complicating the proof. The result given above is sufficient for our purposes, and 
we do not further comment on the possible "optimal" value of the constant p above. Another 
relevant observation is that the inequality in Proposition 3.4 holds true, with the same or 
even smaller p, if we replace Sj^ U Sj^ with a subset of edges S C {Sj^ U Sj^) in (3.16). The 
proof is completely analogous (just U Sj^ is replaced by S). 

4. Preconditioning 

In this section, we present the construction and convergence analysis of the preconditioners 
we propose for the considered IP-methods. 

To construct the preconditioners, we use the subspace splitting given in Proposition 3.1, 
which suggests a simple change of basis. We have that for any u,w & V^^, we can write 
u = z + V, and w = + if, where z,^ & Z and v,ip & V'"^. Therefore, by performing this 
change of basis we can write A{u, w) = A{{z, v), (<^, <j))). The ^o-orthogonality (3.8) of the 
subspaces in the splitting gives 

Ao{{z, v), (C, </))) = Ao{z, + Ao{v, 0). 

which implies that the resulting stiffness matrix of Aq in this new basis is block diagonal. For 
the pure displacement problem (P^r = 0), as discussed in Section 2.4, the spectral equivalence 
given in Lemma 2.2, guarantees that an optimal preconditioner for Aq is also optimal for A. 
Therefore it is enough to study how to efficiently solve each of the blocks in the above block 
diagonal structure of the subproblem resulting from the restriction of Aq to 2 and the 
subproblem on the space V*"^. 

For traction free or mixed type of boundary conditions, although a preconditioner for 
does not result in an optimal solution method. However, the block structure of in the 
new basis already suggests that a reasonable choice for an approximation of ^(■, ■) is 

(4.1) B{{z, «), (C, 0)) = A{z, C) + A{v, cf>). 

The following algorithm describes the application of a preconditioner, which is based on the 
bilinear form in the equation (4.1). 

Algorithm 4.1. Let r e [L'^i^)]'^ be given. Then the action of the preconditioner on r is 
the function u G V^'^ which is obtained from the following three steps. 
1. Find z ^ Z such that 

^(z,C) = (r,C)r, for all C^^Z. 
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2. Find v eV^^ such that 

A{v,cp) = {r,<p)r, for all <^ G F^^. 

3. Set u = z + V . 

As before, the application of this preconditioner corresponds to solving the subproblem of 
the restriction of A{-^ ■) to Z and the subproblem of the restriction of A{-^ ■) to V^^. 

We now briefly discuss how the two smaller sub-problems can be efficiently solved in both 
cases: (1) the case of Dirichlet boundary conditions on all of dVL\ and (2) the case of Neumann 
or mixed boundary conditions. 

Solution in the subspace Z: Lemma 3.3 guarantees that the restriction of A{-, ■) and 
^o('7 ■) to Z is well-conditioned with respect to both, the mesh size and the Lame constants 
A, jjL. Therefore, the linear system corresponding to the subproblem of the restriction to Z 
can be efficiently solved by the method of Conjugate Gradients (CG). A simple consequence 
of the well known estimate on the convergence of CG (see, e.g., [18, 16]) shows that the 
number of CG iterations required to achieve a fixed error tolerance is uniformly bounded, 
independently of the size of the problem and the parameters. 
Solution in V^r. 

We now briefly discuss how to construct a uniform preconditioner for the corresponding 
subproblem on the space V'^^. Rather than developing a completely new method, the idea is 
to use the optimal preconditioners that have already been studied in literature, and modify 
them if needed so that they flt in the present framework. For our discussion, we distinguish 
two cases: the pure displacement problem (Ftv = 0) and the case with mixed or traction free 
boundary conditions {Tn 7^ 0)- 

• For the case of Dirichlet boundary conditions on the entire boundary-the so-called 
pure displacement problem-it is known how to construct optimal order multilevel 
preconditioners that are robust with respect to the parameter A, see e.g. [5, 17, 13] 
and the references therein. 

• The traction free problem or the case of mixed boundary conditions is more difficult 
to handle because the (discrete) Korn inequality is not satisfied for the standard 
discretization by Crouzeix-Raviart elements without additional stabilization, as was 
shown in [11]. The design of optimal and robust solution methods for stabilized 
discretizations is still an open problem, however, auxiliary space techniques might 
bridge this gap soon. 

4.1. Convergence Analysis. We now prove that the proposed block preconditioners are 
indeed optimal so that their convergence is uniform with respect to mesh size and the Lame 
parameters. This result is given in Theorem 4.3. The following Lemma is crucial for this 
proof, since it gives estimates on the norm of the off-diagonal blocks in the 2x2 block 
form of the stiffness matrix associated to ^(•, •), corresponding to the space splitting V^'^ = 
V*"^ © Z. The result provides a measure of the angle between the subspaces and Z, 
with respect to the ^-norm. The proof of this result uses Proposition 3.4. 

Lemma 4.2. Strengthened Cauchy-Schwarz inequality: The following inequality holds 
for any z ^ Z and any v ^ V^^ 

A{z, vf < -i^A{z, z)A(v, v) 

where 7 < 1 and 7 depends only on ao, ai and the constant from Proposition 3.4- 
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Proof. We know that we can always choose large enough, such that for all u G V^*^ we 
have 

Ao{u,u) = {Cs{u) : e{u))r, - 2{i{Ce{u))n}, M)^:.u^:f + «o%,o(M, M) > 0. 

Then it is sufficient to prove that there exists 7 = 7(tti) < 1 such that for aA\ z E Z and for 
all V G the inequality 

holds. By the definition of the spaces 2 and V*"^, on the boundary edges E E Ej^ we have 
either P° [z] = 
we conclude that 



either P° [z] = (if E G or P° |v] = (if ^ G ^^). Hence, from the symmetry of 



[ {M,VUvj)= [ {VUzlM)=0, for all E E Sl and all z E t; G V^^. 
Je J e 

Since for the interior edges E E Sj^ we also have = 0, the above relation and the 

definition of altogether imply that for aA\ z E 2, and v E V'"^ 

(4.2) /(/^e'^^W,H)=«i/3i E /(^^'W,^^M)=0. 

The equation (4.2) and the Schwarz inequality then lead to 



E j ^hflzlM 



-1 2 



-1 2 



< %i(N,H) 



las 



E€£°U£^ 



Next, the result in Proposition 3.4 (more precisely its vector valued form) implies that 



1/2 n 



\0,E 



Ee£°yj£^ 



Therefore, we have 



K,,(W,H)]2 < (i-^)a,,i(H,H) 



/3i E 11^: 



-1/2™ 



\0,E 



< (i--)«.,i(M,W)«.,i(M,M), 



which shows the desired inequality. 



□ 



We are now in a position to prove that the preconditioner given by Algorithm 4.1 is 
uniform with respect to the mesh size and the problem parameters. 
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Theorem 4.3. Let A{-,-) be the symmetric bilinear form defined by (2.24) where 6 = —1 
and ■) be the bilinear form defined by (4-1) ■ Then the following estimates hold for all 
z ^ Z and for all v E V'^^ 

(4.3) j^A{{z, v), (z, v)) < Biiz, v), (z, v)) < -^AHz, v), (z, v)). 

1+7 1-7 

The constant 7 < 1 zs the constant from Lemma 4-2. 

Proof. Using Lemma 4.2 we have 

-2-f^A{z,z)A{v,v) < 2A{z, v) < 2jy/A{z,z)A{v,v) 

and since — — < 2ab < a'^ + b'^ for any real numbers a and b we obtain 

(1 - 7) {A{z, z) + A{v, v)) < A{z, z) + A{v, v) + 2A{z, i;) < (1 + 7) {A{z, z) + A{v, v)) 

which is the same as 

(1 - ^)B{{z, v), (z, v)) < Aiiz, v), (z, v)) < (1 + ^)Biiz, v), (z, v)) 

and thus (4.3) holds with the same constant 7 < 1 as used in the estimate of Lemma 4.2. □ 

Remark 4.4. Note that 7 < g < 1 is uniformly bounded away from 1 and this bound holds 
independently of the parameters h, \, and fi. 

5. Numerical experiments 

In this section we present a set of numerical tests that illustrate our theoretical results. 
We consider the SIPG discretization of the model problem (2.3) on the unit square in IR^ 
with mixed boundary conditions. For the penalty parameters in (2.20) we choose the values 
ao = 4 and ai = 1. The coarsest mesh (at level 0) consists of eight triangles and is refined 
four times. Each refined mesh at level £, £ = 1,2,3,4 is obtained by subdividing every 
triangle at level {£ — 1) into four congruent triangles. The CBS constants and the spectral 
condition numbers summarized in the tables below have been computed using MATLAB. 

In Table 5.1 we list the values of the constant 7^ in the inequality stated in Lemma 4.2 for 
different levels of refinement. Evidently, 7 is uniformly bounded with respect to the mesh 
size (or the number of refinement levels) and also with respect to the material parameters. 
Young's modulus (B and Poisson ratio u (see Remark 4.4). It can be seen from Table 5.2 



Table 5.1. Observed CBS constant 7^ for = (0, 1)^. 



7^ 


V = 0.25 


z/ = 0.4 


u = 0.49 


u = 0.499 




z/ = 0.49999 


£ = 1 


0.0664 


0.025 


0.0024 


2.4024x10^ 


-4 


2.4015x10"*^ 


£ = 2 


0.0678 


0.0255 


0.0025 


2.4567x10" 


-4 


2.4559x10-6 


£ = 3 


0.0684 


0.0258 


0.0025 


2.4866x10" 


-4 


2.4857x10-6 


£^4 


0.0686 


0.0259 


0.0025 


2.4974x10" 


-4 


2.4966x10-6 



that the two subspaces and 2 remain nearly ^-orthogonal when we introduce a jump 
in the Poisson ratio (on the coarsest mesh); In our experiment we set u = Ui = 0.3 (and 
E = El = 1) in the subdomain Qi = [0,0.5] x [0,0.5] U [0.5, 1] x [0.5, 1], and u = U2 (and 
i?2 = 1) in the subdomain Q2 = ^\^2, respectively. 

Next we consider an L-shaped domain Q = [0, 1] x [0, 1] \ (0.5, 1] x (0.5, 1] with Neumann 
boundary conditions on the sides y = and y = I and Dirichlet boundary conditions on 
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the remaining part of the boundary. The initial triangulation (level 0) consists of 4 similar 
triangles. The angle is almost the same as for the square domain, see Table 5.3. 

Furthermore, we computed the relative condition number of the preconditioner B corre- 
sponding to the bilinear form (4.1) for the model problem on the L-shaped domain. The 
results of this experiment, which are listed in Table 5.4, confirm the uniform bound provided 
by Theorem 4.3. 

Finally, we computed the condition number k{Azz) of the matrix A^z related to the re- 
striction of .4(-, ■) to the space Z, again for the model problem on the L-shaped domain. In 
view of Lemma 3.3 we already know that Azz is well-conditioned, and this is clearly seen in 
Table 5.5 where the values of k{Azz) are listed. 



Table 5.2. Observed CBS constant 7^ for Q = (0, 1)^ and jumps in u. 



1' 


z/2 = 0.3 


U2 = 0.4 


U2 = 0.49 


z/2 = 0.499 


z/2 = 0.49999 


i = \ 


0.0451 


0.0177 


0.0442 


0.0509 


0.0517 


1 = 2 


0.0460 


0.0180 


0.0689 


0.0803 


0.0816 


i = 3 


0.0464 


0.0182 


0.0689 


0.0802 


0.0816 


i = 4 


0.0466 


0.0182 


0.0689 


0.0802 


0.0816 



Table 5.3. Observed CBS constant 7^ for L-shaped domain. 



1' 


z/ = 0.25 


u = 0.4 


u = 0.49 


u = 0.499 




z/ = 0.49999 


i = 1 


0.0561 


0.0202 


0.0019 


1.8918x10" 


-4 


1.8906x10-*^ 


i = 2 


0.0631 


0.0233 


0.0022 


2.2118x10- 


-4 


2.2106x10-6 


i = 3 


0.0672 


0.0252 


0.0024 


2.4216x10" 


-4 


2.4207x10-6 


£ = 4 


0.0682 


0.0257 


0.0025 


2.4810x10" 


-4 


2.4801x10-6 



Table 5.4. Tabulated values of k,{B ^A) for L-shaped domain. 



k{B-'A) 


u = 0.25 


u = 0.4 


z/ = 0.49 


u = 0.499 


z/ = 0.49999 


i = l 


1.6204 


1.3314 


1.0912 


1.0279 


1.0028 


i = 2 


1.6713 


1.3606 


1.0990 


1.0302 


1.0030 


i = 3 


1.6997 


1.3774 


1.1037 


1.0316 


1.0031 


i = 4 


1.7073 


1.3820 


1.1050 


1.0320 


1.0032 



Table 5.5. Values of k,{Azz) for L-shaped domain. 



^{Azz) 


V = 0.25 


z/ = 0.4 


u = 0.49 


z/ = 0.499 


jj = 0.49999 


i = l 


8.9067 


7.1484 


6.4788 


6.4220 


6.4158 


i = 2 


9.0875 


7.1932 


6.4829 


6.4229 


6.4164 


i = 3 


9.1577 


7.2080 


6.4841 


6.4230 


6.4164 


e = 4 


9.1794 


7.2118 


6.4844 


6.4230 


6.4164 
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Appendix A. Auxiliary results 

A.l. Bounds on the cardinality of Mi^E) and M2{E). We first recall the definitions of 
No{E), Ui{E) and U2{E), already given in §3.1: 

Mo{E) ■= {T e Th, such that E G T}, 

ATiiE) := {E'eSh, such that 7Vo(E) n 7Vo(^') 7^ 0}> 

M2{E) := {E'eSh, such that ^i{E) f] Mi{E') ^ dl}. 

In the proof of the strengthened Cauchy-Schwarz inequality §3.1 we needed several estimates 
on the cardinality of these sets and these estimates are given in the proposition below. We 
remind the reader that we have \N'o{E)\ < 2. 

Proposition A.l. The following inequalities hold: 

(A.l) \Afi{E)\ < {2d + 1) and \Af2{E)\ < {2d + if . 

Proof. Let E & Shhe fixed. To prove the bound on \Mi{E) \ we consider the elements T G Th, 
such that E & T. In each such element T, there are exactly d faces E' & T, E' ^ E. Since 
there are at most two elements T & Th containing E we have at most 2d faces E' G Sh such 
that E' G Afi{E), and E' 7^ E. Adding E itself to the total count gives \Afi{E)\ <{2d+l). 
The second bound given in (A.l) follows from the first and the following inclusion: 

M2{E)c U Af,{E'). 
E'eJVi{E) 

To show the above inclusion, we consider an arbitrary E" G N'2{E). By the definition of 
J\f2{E), the intersection ofMi{E") and A/i(-B) is not empty. Equivalently, there exists E' G Sh 
such that E' G Afi{E") and E' G Afi{E). On the other hand, from the definition of A/^S"), 
we have that E' G Mi{E") implies that E" G Mi{E'), i.e., if E' is a neighbor of E", then E" 
is a neighbor of E'. 

Putting this together, we conclude that: if E" G M2{E), then there exists E' G Mi{E), 
such that E" G A/i(-E"), and this is exactly the inclusion we wanted to show. 
To prove the desired bound is then straightforward: 



U -^1(^0 

E'eA/'i(E) 



E'eMi{E) E'eJViiE) 
= {2d+l)\Afi{E)\<{2d+l)''. 



□ 
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A.2. A multiplicative relation. This is to prove a basic relation used to derive (3.3) as 
well as (2.15). Let be a map V x W ^ U, where U, V, and W are linear vector spaces 
over the real numbers. We assume that satisfies the following distributive laws: 

aQ{b + c) = aQb + aQc, {a + b)Qc = aQc + bQc, 

and we assume that for all G IR and all ?7 G IR, we have: 

(A.2) i^a)Qirjb) = {^rj){aQb). 

We have the following identities, based on the definitions (2.7): 

(A.3) a+ 6+ - a- r = |a] {{b}} + §a}} {bj. 

Proving this relation is indeed trivial. Some examples for which the reader should verify 
these identities are: (1) For real numbers a and b one may take as the usual multiplication 
of real numbers; (2) a and b elements of a real Hilbert space and inner product; (3) a 
and b are linear operators, and is then the multiplication of linear operators. Note that 
in such case is not necessarily commutative; (4) a is a matrix and 6 is a vector, or more 
generally, a is a linear operator and b is an element of a Hilbert space. 
From (2.7), we have that the right side of the identity (A.3) is 

H {{b}} + {{a}} m = (a+ - a-) (^-^^] + (^^^] {b^ - r ) 



Using the distributive law, and (A.2) (linearity of with respect to scalar multiplication), 
we have 



2 J \ 2 ^ 

_ a-) (6+ + r) + i(a+ + a") (6+ - b~) 

(6+ + r) - ^a- (6+ + r) + (6+ - 6") + ^a" (6+ 

-a+ 6+ + -a+ 6~ - -a" 6+ - -a^ b~ 
2 2 2 2 

1111 
+-a+ 6+ - -a+ 6" + -a" 6+ - -a" 6" 

^ Zi 

1111 

-a+ 6+ - -a" 6~ + -a+ 6+ - -a~ 6~ = a+ 6+ - a" 6". 
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